/* Copyright (c) Colorado School of Mines, 2010.*/
/* All rights reserved.                       */

/* VIEWER3: $Revision: 1.5 $ ; $Date: 1997/07/31 17:10:19 $	*/

#include "par.h" /*this line must appear first*/
#include "GL/glut.h"
#include "MGL/trackball.h"
#include "MGL/viewer3.h"

/*********************** self documentation *****************************/
char *sdoc[] = {
"									",
" VIEWER3 - VIEWER for 3-dimensional model				", 
"									",
" viewer3 [parameters] < hzfile		 				",
" 									",
" Optional Parameters:							",	
" 									",
" hue=1		=1 for hue =0 for black/white				",
" q=-0.6,0.06,-0.06,0.8 define the quaternion			   	",
" tbs=0.8 	the lager the trackball size is, the slower it rotates  ",
" verbose=0     =1 print some useful information			",
" rayfile=NULL	ray path file to read in 				",
" wffile=NULL	wavefront file to read in				",
" sttfile=NULL  surface traveltime file to read in                      ",
"									",
" Menu: The rightmost mouse button, with cursor placed on the graphics  ",
"       window activates the following menu:				",
" Quit:			quits the application				",
" Full Screen:		expands graphics window to full screen		",
" White/Color Rays:		selects colored (default) or white rays	", 
" Plot Rays:			show the rays (default: off)		",
" Surface Traveltimes:	show surface traveltimes (default: off)		",
" Wired or Solid WFs:	style of wavefront display wireframe (default)	",
"			or solid					",	
" Plot Wavefronts:	turn on wavefronts default: off			",
" TRI or TETRA or LAYER or HORZ:	toggles display from wireframe	",
" 			triangulated horizons, to wireframe tetrahedra, ",
"			to solid layers, to solid horizons		",
" Layer1:		toggle layer on or off (default is on)		",
" ....									",
" 									",
" Notes:								",
" In general, the tetrahedra model should be generated by tetramod	",
" (tetrahedra modle builder) which outputs a file called hzfile, so	",
" that viewer3 can read input from stdin.				",
"									",
" If, in addition, you have 3D raypath information (rayfile), wavefronts",
" (wffile) and surface traveltimes (sttfile), typically generated by    ",
" sutetraray (tetrahedral wavefront construction ray tracing),	  	",
" then rays, wavefronts and surface traveltimes can be displayed as well",
" as the model.							 	",
"									",
" The plot may be rotated by depressing the leftmost mouse button	",
" and dragging the cursor in the desired direction of rotation. The speed",
" of rotation is controlled by the speed of the cursor.			",
" 									",
" The plot may be rescaled by depressing the shift key, while also	",
" depressing and dragging the cursor up or down on the display window.	",
" 									",
NULL};

/*
 * Credits:
 *  	CWP: Zhaobo Meng, 1996
 */
/**************** end self doc *******************************************/

int 
main(int argc, char **argv)
{
      int ix;	/*index for nx*/
      int iy;	/*index for ny*/
      int nx;	/*number of sampels in horizon*/
      int ny;	/*number of sampels in horizon*/

      float xmin,xmax;
      float ymin,ymax;
      float zmin,zmax;

      float ***databot;	/*data for plotting*/
      float ***datatop;
      float ***emisbot; /*color on top horizon*/
      float ***emistop; /*color right above base horizon*/
      float v0;

      int verbose;    /*if =1 print some useful information*/

      float eyez;

      int ihz;	/*index for interfaces*/

      int *ntris;	/*number of triangles*/
      int nt;		/*number of samples in each ray*/
      int iray; 	/*index for nrays*/
      int it;		/*index for nt*/
      int iw,iwf,nwf;
      int is,ns;        /*number of sources*/

      float q0[4];

      char *rayfile=""; /*ray file*/
      char *wffile="";
      char *sttfile="";

      FILE *rayfp=NULL;
      FILE *wffp=NULL;
      FILE *sttfp=NULL;

      Layer *horz;
      Layer *ray;
      Layer *wf;
      float vmin=99999.0;
      float vmax=0.0;

      float tt; /*debugging information in the ray file*/

      int itri;

      char names[10];

      int iflag;	/*flag: =1 means ray effective*/
      float emission[4];
      float tmax=0.0,tmin=FLT_MAX;
      float ***stt=NULL;
      float **ttt=NULL;
      int ntr;

      /* hook up getpar */
      initargs(argc,argv);
      requestdoc(1);

      /* get parameters */
      if (!getparint("verbose",&verbose)) verbose=0;

      /******************************************
      Read model parameters from hzfile
      ******************************************/
      fread(&nhz,sizeof(int),1,stdin);
      fread(&nx,sizeof(int),1,stdin);
      fread(&ny,sizeof(int),1,stdin);
      fread(&xmin,sizeof(float),1,stdin);
      fread(&xmax,sizeof(float),1,stdin);
      fread(&ymin,sizeof(float),1,stdin);
      fread(&ymax,sizeof(float),1,stdin);
      fread(&zmin,sizeof(float),1,stdin);
      fread(&zmax,sizeof(float),1,stdin);

      if (verbose)
            fprintf(stderr,"xmin=%f\nxmax=%f\nymin=%f\nymax=%f\nzmin=%f\nzmax=%f\n",
	          xmin,xmax,ymin,ymax,zmin,zmax);

      if (getparstring("rayfile",&rayfile))  
	    if ((rayfp=fopen(rayfile,"r"))==NULL)
		  err("Can not open rayfile %s",rayfile);

      if (getparstring("wffile",&wffile))
	    if ((wffp=fopen(wffile,"r"))==NULL)
		  err("Can not open wffile %s",wffile);

      if (getparstring("sttfile",&sttfile))
	    if ((sttfp=fopen(sttfile,"r"))==NULL)
		  err("Can not open sttfile %s",sttfile);

      if (!getparfloat("tbs",&tbs))   tbs=0.8;
      if (!getparint("hue",&glb_hue)) glb_hue=1; /*1 for glb_hue*/

      if (verbose) 
	    warn("nhz=%d, nx=%d, ny=%d\n",nhz,nx,ny);

      glb_on_or_off=(enum On_or_Off *)ealloc1int(3*nhz+6);
      for (ihz=0;ihz<nhz;ihz++) glb_on_or_off[ihz]=ON;

      horz=(Layer *)alloc1float(sizeof(Layer)*(nhz+1));

      /*********************************************************
      Do not use GLUT_INDEX, which gives no image;
      GLUT_SINGLE will cause redrawing every time you retate it;
      *********************************************************/
      glutInit(&argc, argv);
      glutInitWindowSize(768,768);
      glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH);
      glutCreateWindow("viewer3");
      glutDisplayFunc(redraw);
      glutIdleFunc(NULL);

      if (!getparfloat("q",q0)){
	    q0[0]=-0.6; 
	    q0[1]=0.05;
	    q0[2]=-0.06;
	    q0[3]=0.8;
      }

      normalize_quat(q0);

      curquat[0]=q0[0];
      curquat[1]=q0[1];
      curquat[2]=q0[2];
      curquat[3]=q0[3];

      glutReshapeFunc(myReshape);
      glutVisibilityFunc(vis);
      glutMouseFunc(mouse);
      glutMotionFunc(motion);
      glutCreateMenu(controlLights);
      glutAddMenuEntry("Quit",-1);
      glutAddMenuEntry("Full Screen",0);      
      glutAddMenuEntry("White/Color Rays",1);
      glutAddMenuEntry("Plot Rays",2);
      glutAddMenuEntry("Surface Traveltimes",3);
      glutAddMenuEntry("Wired or Solid WFs",4);
      glutAddMenuEntry("Plot Wavefronts",5);
      glutAddMenuEntry("TRI or TETRA or LAYER or HORZ",6);
      for (ihz=0;ihz<nhz;ihz++) {
	    sprintf(names,"Layer %d",ihz+1);
	    glutAddMenuEntry(names,ihz+7);
      }
      glutAttachMenu(GLUT_RIGHT_BUTTON);

      glShadeModel(GL_SMOOTH);
      glEnable(GL_DEPTH_TEST);
      glEnable(GL_LIGHTING);

      eyez=25;
      glMatrixMode(GL_PROJECTION);
      gluPerspective(
            40.0,   /*fovy: view angle in y direction*/
            1.0,    /*aspect: ratio of width (x) to y (height)*/
            eyez-DIAMETER,  /*near clipping plane*/
            eyez+DIAMETER); /*far clipping plane*/

      glMatrixMode(GL_MODELVIEW);
      gluLookAt(
            0.0, 0.0, eyez, /*(eyex,eyey,eyez): the eye position*/
            0.0, 0.0, 0.0,  /*(centerx,centery,centerz): the center*/
            0.0, 1.0, 0.0); /*(upx,upy,upz): the up direction*/
      glPushMatrix();

      /*the order that tetramod uses is like this*/
      for (ihz=0;ihz<nhz;ihz++) {
            fprintf(stderr,"reading horizon information %d\n",ihz);
	    /**********************************************************
	    input the horizon information from file hzfile:
	    **********************************************************/
	    horz[ihz].x=ealloc2float(nx,ny);
	    horz[ihz].y=ealloc2float(nx,ny);
	    horz[ihz].z=ealloc2float(nx,ny);
	    horz[ihz].v0=ealloc2float(nx,ny);
	    horz[ihz].v1=ealloc2float(nx,ny);

            fprintf(stderr,"read horz[%d].x...\n",ihz);
	    if (fread(horz[ihz].x[0],sizeof(float),nx*ny,
		  stdin)!=nx*ny)
		  err("Can not read x to stdin");
           
            fprintf(stderr,"read horz[%d].y...\n",ihz);
	    if (fread(horz[ihz].y[0],sizeof(float),nx*ny,
		  stdin)!=nx*ny)
		  err("Can not read y to stdin");

	    fprintf(stderr,"read horz[%d].z...\n",ihz);
	    if (fread(horz[ihz].z[0],sizeof(float),nx*ny,
		  stdin)!=nx*ny)
		  err("Can not read z to stdin");

	    fprintf(stderr,"read horz[%d].v0...\n",ihz);
	    if (fread(horz[ihz].v0[0],sizeof(float),nx*ny,stdin)!=
		  nx*ny)
		  err("Can not read v0 to stdin");
 
	    fprintf(stderr,"read horz[%d].v1...\n",ihz);
	    if (fread(horz[ihz].v1[0],sizeof(float),nx*ny,stdin)!=
	      	  nx*ny)
		  err("Can not read v1 to stdin");

	    for (iy=0;iy<ny;iy++) {
	          for (ix=0;ix<nx;ix++) {
		        vmin=MIN(vmin,horz[ihz].v0[iy][ix]);
			vmax=MAX(vmax,horz[ihz].v0[iy][ix]);
		        vmin=MIN(vmin,horz[ihz].v1[iy][ix]);
			vmax=MAX(vmax,horz[ihz].v1[iy][ix]);
		  }
	    }
      }

      if (verbose)
	    fprintf(stderr,"vmin=%f, vmax=%f\n",vmin,vmax);

      horz[nhz].x=ealloc2float(nx,ny);
      horz[nhz].y=ealloc2float(nx,ny);
      horz[nhz].z=ealloc2float(nx,ny);

      fprintf(stderr,"assign horz[%d].x,y,z\n",nhz);
      for (ix=0;ix<nx;ix++) {
	    for (iy=0;iy<ny;iy++) {
	       	  horz[nhz].x[iy][ix]=horz[nhz-1].x[iy][ix];
		  horz[nhz].y[iy][ix]=horz[nhz-1].y[iy][ix];
		  horz[nhz].z[iy][ix]=zmax;
	    }
      }

      databot=ealloc3float(3,nx,ny);
      emisbot=ealloc3float(4,nx,ny);
      datatop=ealloc3float(3,nx,ny);
      emistop=ealloc3float(4,nx,ny);

      for (ihz=0;ihz<nhz;ihz++) {
	    fprintf(stderr,"assigning datatop for ihz=%d\n",ihz);
	    for (ix=0;ix<nx;ix++) {
		  for (iy=0;iy<ny;iy++) {

			datatop[iy][ix][0]=(
			      (horz[ihz].x[iy][ix]-xmin)/
			      (xmax-xmin)-0.5)*DIAMETER;
			datatop[iy][ix][1]=-(
			      (horz[ihz].y[iy][ix]-ymin)/
			      (ymax-ymin)-0.5)*DIAMETER;
			datatop[iy][ix][2]=(
			      (horz[ihz].z[iy][ix]-zmin)/
			      (zmax-zmin)-0.5)*DIAMETER;
			v0=horz[ihz].v0[iy][ix];

			vEmission(v0,vmin,vmax,emistop[iy][ix]);
		  }
	    }

	    fprintf(stderr,"assigning databot for ihz=%d\n",ihz);
	    for (ix=0;ix<nx;ix++) {
		  for (iy=0;iy<ny;iy++) {
			databot[iy][ix][0]=(
			      (horz[ihz+1].x[iy][ix]-xmin)
			      /(xmax-xmin)-0.5)*DIAMETER;
			databot[iy][ix][1]=-(
			      (horz[ihz+1].y[iy][ix]-ymin)
			      /(ymax-ymin)-0.5)*DIAMETER;	
			databot[iy][ix][2]=(
			      (horz[ihz+1].z[iy][ix]-zmin)
			      /(zmax-zmin)-0.5)*DIAMETER;
                                
			v0=horz[ihz].v1[iy][ix];
			vEmission(v0,vmin,vmax,emisbot[iy][ix]);
		  }
	    }

	    showLayer(ihz,databot,datatop,nx,ny,emisbot,emistop);
	    showHorz(ihz,datatop,nx,ny,emistop);
	    showTetra(ihz,databot,datatop,nx,ny,emisbot,emistop);
	    showTri(ihz,datatop,nx,ny,emistop);
      }

      free3float(databot);
      free3float(datatop);
      free3float(emisbot);
      free3float(emistop);

      /*******************************************************************
      The ray positions are generated by sutetraray, named by rayfile. 
      This part will be ignored if rayfile not specified.
      ********************************************************************/
      if (rayfp!=NULL) {

            fscanf(rayfp,
                 "%d =Number of shots\n",&ns);
            fprintf(stderr,"ns=%d\n",ns);

            if (ns<=0 || ns>100) {
                 ns=0;
                 rayfp=NULL;
            }

            ray=(Layer *)alloc1float(sizeof(Layer)*ns);

	    tmax=0.0;
            for (is=0;is<ns;is++) {
                  fscanf(rayfp,
                        "%d =Maximum number of segments\n",&nt);

                  fprintf(stderr,"%d =Maximum number of segments\n",nt);
                  fscanf(rayfp,
                        "%d =Number of rays\n",&ray[is].nrays);
		  fprintf(stderr,"%d =Number of rays\n",ray[is].nrays);

	          ray[is].x=ealloc2float(ray[is].nrays,nt);
	          ray[is].y=ealloc2float(ray[is].nrays,nt);
	          ray[is].z=ealloc2float(ray[is].nrays,nt);
                  ray[is].v0=ealloc2float(ray[is].nrays,nt);
		  ray[is].nseg=ealloc1int(ray[is].nrays);

                  for (iray=0;iray<ray[is].nrays;iray++) {
                        fscanf(rayfp,"%d=nseg %f=ttotal\n",&ray[is].nseg[iray],&tt);

                        if (nt<ray[is].nseg[iray]) err("nt should >=ray[is].nseg[iray]");
                        for (it=0;it<ray[is].nseg[iray];it++) {
                              fscanf(rayfp,"%f %f %f %f %f\n",
                                    &ray[is].x[it][iray],
                                    &ray[is].y[it][iray],
                                    &ray[is].z[it][iray],
				    &ray[is].v0[it][iray],&tt);
			      tmax=MAX(tmax,ray[is].v0[it][iray]);
                        }
                        
                        ray[is].z[ray[is].nseg[iray]-1][iray]=
			      MAX(0.001,ray[is].z[ray[is].nseg[iray]-1][iray]);
           
	                for (it=0;it<ray[is].nseg[iray];it++) {
		      	      ray[is].x[it][iray]=((ray[is].x[it][iray]-xmin)/
		       	            (xmax-xmin)-0.5)*DIAMETER;
		    	      ray[is].y[it][iray]=-((ray[is].y[it][iray]-ymin)/
		       	            (ymax-ymin)-0.5)*DIAMETER;
		       	      ray[is].z[it][iray]=((ray[is].z[it][iray]-zmin)/
		       	            (zmax-zmin)-0.5)*DIAMETER;
	                }
                  }
            }
	    fclose(rayfp);

            /*white rays*/
            glNewList(nhz*4+3,GL_COMPILE);
            emission[0]=emission[1]=emission[2]=emission[3]=1.0;
            glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emission);
            for (is=0;is<ns;is++) {
                  for (iray=0;iray<ray[is].nrays;iray++) {
                        iflag=0;
	                glBegin(GL_LINE_STRIP);
	                for (it=0;it<ray[is].nseg[iray];it++) {
		              if (fabs(ray[is].x[it][iray])<RADIUS &&
		                  fabs(ray[is].y[it][iray])<RADIUS &&
		                  fabs(ray[is].z[it][iray])<RADIUS) {
			            glVertex3f(ray[is].x[it][iray],ray[is].y[it][iray],
			                  ray[is].z[it][iray]);

                                    iflag=1;
		              } else if (iflag) break; /*once good, now bad*/    
	                }
	                glEnd();
                  }
            }
            glEndList();

            /*colored rays*/
            glNewList(nhz*4+4,GL_COMPILE);
            for (is=0;is<ns;is++) {
                  for (iray=0;iray<ray[is].nrays;iray++) {
                        iflag=0;
	                glBegin(GL_LINE_STRIP);
	                for (it=0;it<ray[is].nseg[iray];it++) {
		              if (fabs(ray[is].x[it][iray])<RADIUS &&
		                  fabs(ray[is].y[it][iray])<RADIUS &&
		                  fabs(ray[is].z[it][iray])<RADIUS) {

				    tEmission(
					  ray[is].v0[it][iray],
					  0.0,     /*tmin*/
					  tmax,
					  emission);
				    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emission);
				    glVertex3f(
					  ray[is].x[it][iray],
					  ray[is].y[it][iray],
			                  ray[is].z[it][iray]);
                                    iflag=1;
		              } else if (iflag) break; /*once good, now bad*/    
	                }
	                glEnd();
                  }
            }
            glEndList();
      }

      /*************************************************************
      Plot the wavefront if it is given. If the wffile does not 
      contain effective data, ntris may be wild. In this case, do
      thing about the wavefront.
      *************************************************************/
      if (wffp!=NULL) {
            fscanf(wffp,"%d = nwf2dump\n",&nwf);
            fprintf(stderr,"nwf2dump=%d\n",nwf);
            if (nwf>200) wffp=NULL;
      }

      if (wffp!=NULL) {

   	    emission[0]=1.0;
            emission[1]=1.0; 
	    emission[2]=0.0;
	    emission[3]=1.0;

	    wf=(Layer *)alloc1float(sizeof(Layer)*nwf);
            ntris=ealloc1int(sizeof(int)*nwf);

	    for (iwf=0;iwf<nwf;iwf++) {

	          if (1!=fscanf(wffp,"%d = ntris\n",&ntris[iwf])) {
		        nwf=iwf;
                        break;
	          }

		  if (ntris[iwf]==0) {
                        nwf=iwf;
                        break;
	          }

		  if (verbose)
			warn("ntris=%d of nwf=%d\n",ntris[iwf],nwf);

		  wf[iwf].x=ealloc2float(3,ntris[iwf]);
		  wf[iwf].y=ealloc2float(3,ntris[iwf]);
		  wf[iwf].z=ealloc2float(3,ntris[iwf]);

		  for (it=0;it<ntris[iwf];it++) {
                        fscanf(wffp,"%f %f %f %f %f %f %f %f %f\n",
                              wf[iwf].x[it],  wf[iwf].y[it],  wf[iwf].z[it],
                              wf[iwf].x[it]+1,wf[iwf].y[it]+1,wf[iwf].z[it]+1,
                              wf[iwf].x[it]+2,wf[iwf].y[it]+2,wf[iwf].z[it]+2);
		  }

		  fprintf(stderr,"Totally read in %d wavefront triangles\n",ntris[iwf]);

		  for (it=0;it<ntris[iwf];it++) {
			for (iw=0;iw<3;iw++) {
			      wf[iwf].x[it][iw]=((wf[iwf].x[it][iw]-xmin)/
					(xmax-xmin)-0.5)*DIAMETER;
			      wf[iwf].y[it][iw]=-((wf[iwf].y[it][iw]-ymin)/
					(ymax-ymin)-0.5)*DIAMETER;
			      wf[iwf].z[it][iw]=((wf[iwf].z[it][iw]-zmin)/
					(zmax-zmin)-0.5)*DIAMETER;
			}
       	          }
	    }
	    fclose(wffp);

            fprintf(stderr,"Click right MB to get menu\n");
            fprintf(stderr,"Click left MB and drag to rotate\n");
            fprintf(stderr,"Press shift and push left MB to scale\n");

  	    glNewList(nhz*4+6,GL_COMPILE);
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emission);

	    for (iwf=0;iwf<nwf;iwf++) {
	          for (itri=0;itri<ntris[iwf];itri++) {
		        glBegin(GL_LINE_LOOP);
		        if (fabs(wf[iwf].x[itri][0])<RADIUS &&
		            fabs(wf[iwf].y[itri][0])<RADIUS &&
		            fabs(wf[iwf].z[itri][0])<RADIUS && 
		            fabs(wf[iwf].x[itri][1])<RADIUS &&
		            fabs(wf[iwf].y[itri][1])<RADIUS &&
		            fabs(wf[iwf].z[itri][1])<RADIUS &&
                            fabs(wf[iwf].x[itri][2])<RADIUS &&
                            fabs(wf[iwf].y[itri][2])<RADIUS &&
                            fabs(wf[iwf].z[itri][2])<RADIUS) {
			          glVertex3f(wf[iwf].x[itri][0],
					     wf[iwf].y[itri][0],
					     wf[iwf].z[itri][0]);
			          glVertex3f(wf[iwf].x[itri][1],
					     wf[iwf].y[itri][1],
					     wf[iwf].z[itri][1]);
			          glVertex3f(wf[iwf].x[itri][2],
					     wf[iwf].y[itri][2],
					     wf[iwf].z[itri][2]);
		        } else {
                            fprintf(stderr,"warning: some triangles ignored\n");
                            glEnd();
                            break;
                        }                        
		        glEnd();
	          }
            }
            glEndList();

            /*solid wavefronts*/
  	    glNewList(nhz*4+7,GL_COMPILE);
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emission);

	    for (iwf=0;iwf<nwf;iwf++) {
	          for (itri=0;itri<ntris[iwf];itri++) {
		        glBegin(GL_TRIANGLE_STRIP);
		        if (fabs(wf[iwf].x[itri][0])<RADIUS &&
		            fabs(wf[iwf].y[itri][0])<RADIUS &&
		            fabs(wf[iwf].z[itri][0])<RADIUS && 
		            fabs(wf[iwf].x[itri][1])<RADIUS &&
		            fabs(wf[iwf].y[itri][1])<RADIUS &&
		            fabs(wf[iwf].z[itri][1])<RADIUS &&
                            fabs(wf[iwf].x[itri][2])<RADIUS &&
                            fabs(wf[iwf].y[itri][2])<RADIUS &&
                            fabs(wf[iwf].z[itri][2])<RADIUS) {
			          glVertex3f(wf[iwf].x[itri][0],
					     wf[iwf].y[itri][0],
					     wf[iwf].z[itri][0]);
			          glVertex3f(wf[iwf].x[itri][1],
					     wf[iwf].y[itri][1],
					     wf[iwf].z[itri][1]);
			          glVertex3f(wf[iwf].x[itri][2],
					     wf[iwf].y[itri][2],
					     wf[iwf].z[itri][2]);
		        } else {
                            fprintf(stderr,"warning: some triangles ignored\n");
                            glEnd();
                            break;
                        }                        
		        glEnd();
	          }
            }
            glEndList();
      }

      /*surface traveltimes*/
      if (sttfp!=NULL) {
            fscanf(sttfp,"%d = ntris\n",&ntr);
            fprintf(stderr,"ntr=%d\n",ntr);
            if (ntr>2000) sttfp=NULL;
      }      

      if (sttfp!=NULL && ntr>0) {
            stt=ealloc3float(3,3,ntr);
            ttt=ealloc2float(3,ntr);

            tmax=0.0;
            tmin=1.0e+10;
	    for (itri=0;itri<ntr;itri++) {
	          fscanf(sttfp,"%f %f %f %f %f %f %f %f %f %f %f %f\n",
		       &stt[itri][0][0],
		       &stt[itri][0][1],
		       &stt[itri][0][2],
		       &ttt[itri][0],
		       &stt[itri][1][0],
		       &stt[itri][1][1],
		       &stt[itri][1][2],
		       &ttt[itri][1],
		       &stt[itri][2][0],
		       &stt[itri][2][1],
		       &stt[itri][2][2],
		       &ttt[itri][2]);

                  tmax=MAX(tmax,ttt[itri][0]);
		  tmax=MAX(tmax,ttt[itri][1]);
		  tmax=MAX(tmax,ttt[itri][2]);

                  tmin=MIN(tmin,ttt[itri][0]);
		  tmin=MIN(tmin,ttt[itri][1]);
		  tmin=MIN(tmin,ttt[itri][2]);

		  stt[itri][0][0]=((stt[itri][0][0]-xmin)/
		       (xmax-xmin)-0.5)*DIAMETER;
		  stt[itri][0][1]=-((stt[itri][0][1]-ymin)/
		       (ymax-ymin)-0.5)*DIAMETER;
		  stt[itri][0][2]=((stt[itri][0][2]-zmin)/
		       (zmax-zmin)-0.5)*DIAMETER;

		  stt[itri][1][0]=((stt[itri][1][0]-xmin)/
		       (xmax-xmin)-0.5)*DIAMETER;
		  stt[itri][1][1]=-((stt[itri][1][1]-ymin)/
		       (ymax-ymin)-0.5)*DIAMETER;
		  stt[itri][1][2]=((stt[itri][1][2]-zmin)/
		       (zmax-zmin)-0.5)*DIAMETER;

		  stt[itri][2][0]=((stt[itri][2][0]-xmin)/
		       (xmax-xmin)-0.5)*DIAMETER;
		  stt[itri][2][1]=-((stt[itri][2][1]-ymin)/
		       (ymax-ymin)-0.5)*DIAMETER;
		  stt[itri][2][2]=((stt[itri][2][2]-zmin)/
		       (zmax-zmin)-0.5)*DIAMETER;
            }
      }

      tmax=MAX(tmax,tmin+0.01);

      glNewList(nhz*4+5,GL_COMPILE);
      for (itri=0;itri<ntr;itri++) {
	    glBegin(GL_TRIANGLE_STRIP);

            tEmission(ttt[itri][0],tmin,tmax,emission);
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emission);
	    glVertex3fv(stt[itri][0]);

	    tEmission(ttt[itri][1],tmin,tmax,emission);
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emission);
	    glVertex3fv(stt[itri][1]);

	    tEmission(ttt[itri][2],tmin,tmax,emission);
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emission);
	    glVertex3fv(stt[itri][2]);

	    glEnd();
      }
      glEndList();

      glutMainLoop();
      return 0;
}

/*********************************************************************
Author: Zhaobo Meng, CWP 1996
*********************************************************************/
void
showHorz(

      int ihz,float ***data,
      int nx,int ny,float ***emis)
{
      int ix,iy;
      int nxhalf;

      nxhalf=nx/2;

      glNewList(ihz+nhz+3,GL_COMPILE);
      glBegin(GL_QUAD_STRIP);
      for (ix=0;ix<nx-1;ix++) {
      	    for (iy=0;iy<ny;iy++) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix]);
		  glVertex3fv(data[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix+1]);
		  glVertex3fv(data[iy][ix+1]);
	    }

	    for (iy=ny-1;iy>=0;iy--) {
	       	  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix+1]);
		  glVertex3fv(data[iy][ix+1]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix]);
		  glVertex3fv(data[iy][ix]);
	    }
      }

      glEnd();
      glEndList();
}

void
showLayer(int ihz,float ***databot,float ***datatop,
      int nx,int ny,float ***emisbot,float ***emistop) 
{
      int ix,iy;

      glNewList(ihz+3,GL_COMPILE);
      glBegin(GL_QUAD_STRIP);
      for (ix=0;ix<nx-1;ix++) {
	    for (iy=0;iy<ny;iy++) {
		  if (databot[iy][ix][2]<EPS+datatop[iy][ix][2] &&
		      databot[iy][ix+1][2]<EPS+datatop[iy][ix+1][2]) continue; 
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix]);
		  glVertex3fv(databot[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix+1]);
		  glVertex3fv(databot[iy][ix+1]);
	    }
	    for (iy=ny-1;iy>=0;iy--) {
		  if (databot[iy][ix][2]<EPS+datatop[iy][ix][2] &&
		      databot[iy][ix+1][2]<EPS+datatop[iy][ix+1][2]) continue; 
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix]);
		  glVertex3fv(datatop[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix+1]);
		  glVertex3fv(datatop[iy][ix+1]);
	    }
      }

      for (iy=0;iy<ny;iy++) {
	    if (databot[iy][nx-1][2]<EPS+datatop[iy][nx-1][2]) continue;
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][nx-1]);
	    glVertex3fv(databot[iy][nx-1]);
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][nx-1]);
	    glVertex3fv(datatop[iy][nx-1]);
      }
      for (ix=nx-1;ix>=0;ix--) {
	    if (databot[ny-2][ix][2]<EPS+datatop[ny-2][ix][2]) continue;
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[ny-1][ix]);
	    glVertex3fv(databot[ny-1][ix]);
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[ny-1][ix]);
	    glVertex3fv(datatop[ny-1][ix]);
      }

      for (iy=ny-1;iy>=0;iy--) {
	    if (databot[iy][0][2]<EPS+datatop[iy][0][2]) continue;
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][0]);
	    glVertex3fv(databot[iy][0]);
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][0]);
	    glVertex3fv(datatop[iy][0]);
      }

      for (ix=0;ix<nx;ix++) {
	    if (databot[0][ix][2]<EPS+datatop[0][ix][2]) continue;
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[0][ix]);
	    glVertex3fv(databot[0][ix]);
	    glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[0][ix]);
	    glVertex3fv(datatop[0][ix]);
      }
      glEnd();
      glEndList();
}

void    /*notice: I changed the order of 0 and 1 here to make it upside down*/
showTetra(int ihz,float ***datatop,float ***databot,
      int nx,int ny,float ***emistop,float ***emisbot) 
{
      int ix,iy;

      glNewList(ihz+3+nhz*2,GL_COMPILE);
      glBegin(GL_LINES);
      /**********************************************************
      Plot 2 grids like this

         (ix,iy+2)_____(ix+1,iy+2)
                 |    /|
                 |   / |
	         |  /  |
                 | /   |
                 |/    |
 	(ix,iy+1)|-----|(ix+1,iy+1)
		 |\    |
		 | \   |
		 |  \  |
		 |   \ |
	  (ix,iy)|____\|(ix+1,iy)

      **********************************************************/
      for (ix=0;ix<nx-1;ix++) {   /*horizontal edges*/
	    for (iy=0;iy<ny;iy++) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix]);
		  glVertex3fv(databot[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix+1]);
		  glVertex3fv(databot[iy][ix+1]);
            }
      }

      for (ix=0;ix<nx;ix++) {   /*vertical edges*/
	    for (iy=0;iy<ny-1;iy++) {
	          glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix]);
		  glVertex3fv(databot[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy+1][ix]);
		  glVertex3fv(databot[iy+1][ix]);
            }
      }

      for (ix=0;ix<nx-1;ix++) {   /*slant edges*/
	    for (iy=0;iy<ny-1;iy+=2) {
	          glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy+1][ix]);
		  glVertex3fv(databot[iy+1][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix+1]);
		  glVertex3fv(databot[iy][ix+1]);
	    }
        
	    for (iy=1;iy<ny-1;iy+=2) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix]);
		  glVertex3fv(databot[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy+1][ix]);
		  glVertex3fv(databot[iy+1][ix+1]);
	    }
      }

      /*the base horizon*/
      for (ix=0;ix<nx-1;ix++) {   /*horizontal edges*/
	    for (iy=0;iy<ny;iy++) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix]);
		  glVertex3fv(datatop[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix+1]);
		  glVertex3fv(datatop[iy][ix+1]);
            }
      }

      for (ix=0;ix<nx;ix++) {   /*vertical edges*/
	    for (iy=0;iy<ny-1;iy++) {
	          glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix]);
		  glVertex3fv(datatop[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy+1][ix]);
		  glVertex3fv(datatop[iy+1][ix]);
            }
      }

      for (ix=0;ix<nx-1;ix++) {   /*slant edges*/
	    for (iy=0;iy<ny-1;iy+=2) {
	          glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy+1][ix]);
		  glVertex3fv(datatop[iy+1][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix+1]);
		  glVertex3fv(datatop[iy][ix+1]);
	    }
        
	    for (iy=1;iy<ny-1;iy+=2) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix]);
		  glVertex3fv(datatop[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy+1][ix+1]);
		  glVertex3fv(datatop[iy+1][ix+1]);
	    }
      }

      /**********************************************************
      Plot vertical segments
      **********************************************************/
      for (ix=0;ix<nx;ix++) {
	    for (iy=0;iy<ny;iy++) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix]);
		  glVertex3fv(databot[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix]);
		  glVertex3fv(datatop[iy][ix]);
            }
      }

      /**********************************************************
      Plot slant segments
      **********************************************************/
      for (ix=0;ix<nx-1;ix++) {
	    for (iy=0;iy<ny-1;iy+=2) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix]);
		  glVertex3fv(databot[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix+1]);
		  glVertex3fv(datatop[iy][ix+1]);

		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy+1][ix]);
		  glVertex3fv(databot[iy+1][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix]);
		  glVertex3fv(datatop[iy][ix]);

		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy+1][ix]);
		  glVertex3fv(databot[iy+1][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix+1]);
		  glVertex3fv(datatop[iy][ix+1]);

		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy+1][ix+1]);
		  glVertex3fv(databot[iy+1][ix+1]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy+1][ix]);
		  glVertex3fv(datatop[iy+1][ix]);

		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy+1][ix+1]);
		  glVertex3fv(databot[iy+1][ix+1]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy][ix+1]);
		  glVertex3fv(datatop[iy][ix+1]);
	    }

	    for (iy=1;iy<ny-1;iy+=2) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix]);
		  glVertex3fv(databot[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy+1][ix]);
		  glVertex3fv(datatop[iy+1][ix]);

		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix]);
		  glVertex3fv(databot[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy+1][ix+1]);
		  glVertex3fv(datatop[iy+1][ix+1]);

		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy+1][ix]);
		  glVertex3fv(databot[iy+1][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy+1][ix+1]);
		  glVertex3fv(datatop[iy+1][ix+1]);

		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix]);
		  glVertex3fv(databot[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy+1][ix+1]);
		  glVertex3fv(datatop[iy+1][ix+1]);

		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emisbot[iy][ix+1]);
		  glVertex3fv(databot[iy][ix+1]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emistop[iy+1][ix+1]);
		  glVertex3fv(datatop[iy+1][ix+1]);
	    }
      }

      glEnd();
      glEndList();
}

void
showTri(int ihz,float ***data,
      int nx,int ny,float ***emis) 
{
      int ix,iy;

      glNewList(ihz+3+nhz*3,GL_COMPILE);
      glBegin(GL_LINES);

      /**********************************************************
	Plot 2 grids like this

         (ix,iy+2)_____(ix+1,iy+2)
                 |    /|
                 |   / |
	         |  /  |
                 | /   |
                 |/    |
 	(ix,iy+1)|-----|(ix+1,iy+1)
		 |\    |
		 | \   |
		 |  \  |
		 |   \ |
	  (ix,iy)|____\|(ix+1,iy)

      **********************************************************/
      for (ix=0;ix<nx-1;ix++) {   /*horizontal edges*/
	    for (iy=0;iy<ny;iy++) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix]);
		  glVertex3fv(data[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix+1]);
		  glVertex3fv(data[iy][ix+1]);
            }
      }

      for (ix=0;ix<nx;ix++) {   /*vertical edges*/
	    for (iy=0;iy<ny-1;iy++) {
	          glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix]);
		  glVertex3fv(data[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy+1][ix]);
		  glVertex3fv(data[iy+1][ix]);
            }
      }

      for (ix=0;ix<nx-1;ix++) {   /*slant edges*/
	    for (iy=0;iy<ny-1;iy+=2) {
	          glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix]);
		  glVertex3fv(data[iy+1][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix+1]);
		  glVertex3fv(data[iy][ix+1]);
	    }
        
	    for (iy=1;iy<ny-1;iy+=2) {
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy][ix]);
		  glVertex3fv(data[iy][ix]);
		  glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emis[iy+1][ix+1]);
		  glVertex3fv(data[iy+1][ix+1]);
	    }
      }

      glEnd();
      glEndList();
}

void
recalcModelView(void)
{
      GLfloat m[4][4];		/*the matrix*/

      glPopMatrix();
      glPushMatrix();

      build_rotmatrix(m, curquat);
      glMultMatrixf(&m[0][0]);
      glScalef(glb_scalefactor, glb_scalefactor, glb_scalefactor);
      glb_newmodel=0;
}

/***************************************************************
function to print a message in the window
***************************************************************/
void
showMessage(GLfloat x, GLfloat y, GLfloat z, char *message)
{
      glPushMatrix();
      glDisable(GL_LIGHTING);
      glTranslatef(x, y, z);
      glScalef(.01, .01, .01);
      /* To rotate 90 degrees counterclockwise aroubd (0,0,1)*/
      glRotatef(90.0, 0.0, 0.0, 1.0);
      while (*message) {
       	    glutStrokeCharacter(GLUT_STROKE_ROMAN, *message);
	    message++;
      }
      glEnable(GL_LIGHTING);
      glPopMatrix();
}

/***************************************************
Called whenever glutPostRedisplay() is called
***************************************************/
void
redraw(void)
{
      int ihz;
      if (glb_newmodel) recalcModelView();
      glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

      for (ihz=0;ihz<nhz;ihz++) {
	    if (hlt==LAYER && glb_on_or_off[ihz]==ON)
		  glCallList(ihz+3);
	    else if (hlt==HORZ && glb_on_or_off[ihz]==ON)
		  glCallList(nhz+ihz+3);
	    else if (hlt==TETRA && glb_on_or_off[ihz]==ON)
		  glCallList(ihz+2*nhz+3);
	    else if (hlt==TRI && glb_on_or_off[ihz]==ON)
		  glCallList(ihz+3*nhz+3);
      }

      if (glb_plot_rays==PLOT_RAYS) {
	    if (wc_rays==WHITE)
	          glCallList(nhz*4+3);
	    else
	          glCallList(nhz*4+4);
      }

      if (plot_st==PLOT_SURFACE_TRAVELTIMES)
            glCallList(nhz*4+5);

      if (glb_plot_wf==PLOT_WF) {
	    if (wfs==WIRED)
                  glCallList(nhz*4+6);
	    else
	          glCallList(nhz*4+7);
      }

      glFlush();
      /*showMessage(0, 0, 0, "x direction");*/
      glutSwapBuffers();
}

/****************************************
This call back is called when the window
is rescaled
****************************************/
void
myReshape(int w, int h)
{
      glViewport(0, 0, w, h);
      glb_W=w;
      glb_H=h;
}

/*****************************************
This call back is called when the mouse is
pressed
*****************************************/
void
mouse(int button, int state, int x, int y)
{
      if (button==GLUT_LEFT_BUTTON && state==GLUT_DOWN) {
	    glb_spinning=0;
	    glutIdleFunc(NULL);
	    glb_moving=1;
	    glb_beginx=x;
	    glb_beginy=y;
	    if(glutGetModifiers() & GLUT_ACTIVE_SHIFT) 
		  glb_scaling=1;
	    else 
		  glb_scaling=0;
      }	
      if (button==GLUT_LEFT_BUTTON && state==GLUT_UP) {
	    glb_moving=0;
	    glb_spinning=0;
      }
}

/********************************************
Animate the object to continue its motion
********************************************/
void
animate(void)
{
      add_quats(lastquat, curquat, curquat);
      glb_newmodel=1;
      glutPostRedisplay();
}

/*********************************************
Called when the mouse is pressed while moved
*********************************************/
void
motion(int x,int y)
{
      if (glb_scaling) {
	    glb_scalefactor=glb_scalefactor*
		  (1.0+(((float)(glb_beginy-y))/glb_H));
	    glb_beginx=x;
	    glb_beginy=y;
	    glb_newmodel=1;
	    glutPostRedisplay();
	    return;
      }

      if (glb_moving) {
	    trackball(lastquat,
		  (2.0*glb_beginx-glb_W)/glb_W,
		  (glb_H-2.0*glb_beginy)/glb_H,
		  (2.0*x-glb_W)/glb_W,
		  (glb_H-2.0*y)/glb_H,tbs);

	    #ifdef DEBUG
	    fprintf(stderr,"q=%f,%f,%f,%f\n",curquat[0],curquat[1],
		  curquat[2],curquat[3]);
            #endif

	    glb_beginx=x;
	    glb_beginy=y;
	    glb_spinning=1;
	    glutIdleFunc(animate);
      }
}

/******************************************
switches
******************************************/
void
controlLights(int value)
{
      switch (value) {
	    case -1:
		  exit(0);
		  break;
	    case 0:
		  glutFullScreen();
		  break;
            case 1:
                  if (wc_rays==WHITE)
                        wc_rays=COLORED;
                  else
                        wc_rays=WHITE;
                  break;
            case 2:
		  if (glb_plot_rays==DO_NOT_PLOT_RAYS)
		       	glb_plot_rays=PLOT_RAYS;
		  else
			glb_plot_rays=DO_NOT_PLOT_RAYS;
                  break;
            case 3:
                  if (plot_st==DO_NOT_PLOT_SURFACE_TRAVELTIMES)
		        plot_st=PLOT_SURFACE_TRAVELTIMES;
                  else 
		        plot_st=DO_NOT_PLOT_SURFACE_TRAVELTIMES;
		  break;
            case 4:
                  if (wfs==WIRED)
                        wfs=SOLID;
                  else
	      	        wfs=WIRED;                  
                  break;
            case 5:
                  if (glb_plot_wf==DO_NOT_PLOT_WF)
                        glb_plot_wf=PLOT_WF;
                  else
                        glb_plot_wf=DO_NOT_PLOT_WF;
                  break;
            case 6:
		  if (hlt==HORZ)
		        hlt=TRI;
		  else if (hlt==TRI)
			hlt=TETRA;
		  else if (hlt==TETRA)
			hlt=LAYER;
		  else 
			hlt=HORZ;	
		  break;
	    default:
		  if (glb_on_or_off[value-7]==ON)
			glb_on_or_off[value-7]=OFF;
		  else 
		        glb_on_or_off[value-7]=ON;
	          break;	
      }
      glutPostRedisplay();
}

/************************************************
When it is not visible, do not waste the computer
resources
************************************************/
void
vis(int visible)
{
      if (visible==GLUT_VISIBLE) {
	    if (glb_spinning)
		  glutIdleFunc(animate);
      } else {
	    if (glb_spinning)
	       	  glutIdleFunc(NULL);
      }
}

/*********************************************************************
Author: Zhaobo Meng, CWP 1996

Set up the traveltime mapping: traveltime->color
************************************************/
void
vEmission(float v,float vmin,float vmax,
      float *emission)
{
      float color;
      color=(v-vmin)/(vmax-vmin);
      color=MAX(0.0,MIN(1.0,color));
      if (glb_hue==1) {
	    emission[0]=1.0;
	    emission[1]=color;
	    emission[2]=1.0;
      } else {
	    emission[0]=color;
	    emission[1]=color;
	    emission[2]=color;
      }
      emission[3]=glb_alpha;
}

/*********************************************************************
Author: Zhaobo Meng, CWP 1996

Set up the traveltime mapping: traveltime->color
*********************************************************************/
void
tEmission(float t,float tmin,float tmax,
      float *emission)
{
      if (fabs(tmax-tmin)==0.0)
            emission[1]=0.0;
      else  
            emission[1]=MAX(0.0,MIN(1.0,(t-tmin)/(tmax-tmin)));
      emission[0]=1.0;
      emission[2]=emission[1];
      emission[3]=glb_alpha;
}

static void 
normalize_quat(float q[4])
{
      int i;
      float mag;

      mag = (q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
      for (i = 0; i < 4; i++) q[i] /= mag;
}
